Variational approach to vortex penetration and vortex interaction 
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A variational calculation for vortex penetration is presented. Variational trial functions for the 
Meissner state are combined with variational functions for a vortex near the surface. The latter 
is based on Clem's trial solutions for a vortex in bulk, which were adapted to include surface 
effects through consideration of an image vortex. Three variational parameters are considered, 
corresponding to the effective coherence length of the vortex, the effective penetration length for the 
Meissner currents, and the value of the order parameter at the surface. The results show that the last 
two variational parameters are independent of vortex position. Explicit calculations are presented 
for several k values. The energy barrier for vortex penetration is shown to be in good agreement 
with full numerical calculations of the Ginzburg-Landau equations. We consider the variation of the 
magnetic flux carried by a vortex as it gets inside the superconductor and agreement with known 
experimental and theoretical results is obtained. The model was extended to calculate the force 
between two vortices, and the results show that the force goes to zero as the pair comes close to 
the surface. This result can be of interest for the study of the melting of the vortex lattice and for 
vortices confined in mesoscopic superconductors. The variational approach can be very helpful for 
intermediate k values when numerical calculations become computationally demanding because it 
provides manageable expressions for all physically relevant quantities. 

PACS numbers: 74.20.De, 74.25. Ha, 74.25.Op 



I. INTRODUCTION 

The behavior of a vortex near the surface of a su- 
perconductor has been the subject of several recent 
papers Semi analytical results have been known 
since the Bean-Livingston 7 - model was formulated and 
simple calculations were made from the London model 
by de Gennes^ Among other properties, these calcu- 
lations give the characteristics of the surface barrier 
for vortex penetration. The geometrical surface barrier 
in superconducting thin films has been considered for 
high-T c superconductors in Refs. [l] and [3- While the 
Bean Livingston barrier is of energetic origin, the geo- 
metrical barrier is strongly dependent on sample shape. 
The surface barriers are also very important in meso- 
scopic superconductors'^^ Interesting results pertain- 
ing to the ac magnetic properties of mesoscopic supercon- 
ductors have been obtained from numerical calculations 
based on the finite-difference method in Refs. Hand©. 

On the other hand, variational calculations are known 
to provide manageable and accurate results in many dif- 
ferent physical problems. In other applications of the 
Ginzburg-Landau (GL) equations, variational calcula- 
tions are known to give good agreement with exact re- 
sults. In some cases, variational calculations have pre- 
ceded exact or numerical calculations. Such is the case 
for surface superconductivity^^ the mixed state in type 
II superconductors^ and, more recently, superconduct- 
ing micronet works ^2 

In this paper, we present a variational approach to 
the solution of the GL equations for a vortex near the 
surface of a superconductor, starting from Clem's vari- 
ational ansatsi^ for a vortex in bulk. In this form, we 
are able to compare the variational results to the full nu- 



merical calculations. Vortices appear in the presence of 
an externally applied field, which also induces Meissner 
currents. For this reason, it is necessary to variationally 
model both aspects of the behavior of a superconduc- 
tor. The Clem ansatz has also been used recently in the 
context of mesoscopic superconductors in Ref. Il4 . 

The paper is organized as follows: In Sec. II. A, we 
present a variational description of the Meissner state, 
and the variational solution is compared to the full nu- 
merical results of the GL equations in Sec. II B. In Sec. 
Ill, Clem's ansatz for a single vortex in a bulk material is 
adapted to describe a vortex near the surface of a super- 
conductor and is combined with the description of the 
Meissner state. In Sees. Ill A and III B, we present 
the results of the variational calculations including these 
three parameters: the penetration length for the Meiss- 
ner currents, the order parameter at the sample surface, 
and the coherence length for the vortex size. It turns out 
that the first two parameters are independent of the vor- 
tex position. The results of the variational calculation 
are compared to the full numerical results, particularly 
for the energy barrier for k = 2 and k = 3. Quite rea- 
sonable agreement between both methods is obtained for 
both k values. In Sec. IV, the model is extended to 
calculate the force between two vortices as a function of 
their distance to the surface. We show that the interac- 
tion force goes to zero as the pair approaches the sample 
surface. Finally, in Sec. V, we give our conclusions. 



II. DESCRIPTION OF THE MEISSNER STATE 

In this section, we propose a variational model to de- 
scribe the Meissner state of a semi-infinite sample. In 
Sec. II A we obtain an approximate solution of the GL 
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equations valid at low magnetic fields when depletion of 
the order parameter at the sample surface is small. By 
using this approximate solution, we propose a variational 
model to describe the Meissner state at higher values of 
the field. In Sec. II B, the variational solution is com- 
pared to the full numerical results of the GL equations. 

A. Variational model for the Meissner state 

By writing the order parameter as ip = fe llp , we ob- 
tain the following expression in normalized units for the 
difference between the free energy of the normal and su- 
perconducting states (Q s — Q n ): 



h% = /(/ 2 -i) + (A,)V, (4) 

t£ - > 2A - < 5) 

with the following boundary conditions: (df / dx)\ x= o = 
and B|a; = o = {dA y /dx)\ x=Q = H a . 

An approximate solution of these equations at low 
fields can be found by assuming f(x) = foo — f]( x )i with 
\r]{x)\ -C foo- In the present normalization, f^ = 1; thus, 
we can write the solution to Eq. ([5]) when r)(x) — > as 
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Lengths r are scaled in units of the zero temperature 
penetration length A(0), the externally applied magnetic 
field H a and B also in units of ^/2H C (0), the vector po- 
tential A in units of y/2X(0)H c (Q), the order parameter 
ip in units of -0oo, the current J in units of i&^eh/m^, 
and velocities in units of ft/2m£(0). 

The GL equations become 
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We assume a semi-infinite medium subjected to a mag- 
netic field parallel to the superconductor-vacuum inter- 
face. We choose the x axis perpendicular to this inter- 
face and take the z direction parallel to the applied field 
B = B z (x)z. 

Equations (f5]) and © must be complemented with the 
appropriate boundary conditions at the sample surface, 
which when separated into its real and imaginary parts, 
imply (V/)j_| s = and (u)j_| s = 0, where the first re- 
lation indicates that the slope of the order parameter 
perpendicular to the surface must be zero at the surface, 
whereas the second implies that the velocity of the super- 
conducting electrons (u ~ (V</J — kA)) has no component 
perpendicular to the surface. For a semi-infinite sample 
with no demagnetizing effects, the condition B| s = H a , 
where H a is the externally applied field, also applies at 
the surface. 

In this configuration, the order parameter depends 
only on x, f = f(x). Moreover, V x B has only a 
non- vanishing component, (V x B) = d z B x — d x B z — 
—d x B z (x)y. From Eq. (J3]), we then have A = A y (x)y. 

In the London gauge, the order parameter is real and 
we can eliminate the phase ip in the GL equations. In this 
geometry and specific gauge, Eqs. @ and become 



Ay ~ B a B 



(0) 



In the following, we assume that at low fields for f(x) < 
1, the vector potential can be conveniently approximated 
by a variational expression of the following form: 



A y — —\MH a e 



—x/Xn 



(7) 



where Aj\./ is a variational parameter. As will be the case 
for the other two variational parameters to be introduced 
later, Am is a field and temperature dependent parame- 
ter. 

By using Eq. ([7]) and f(x) = 1 — ij(x), in a first ap- 
proximation, Eq. ([?]) becomes 



1 d 2 T\ _ 2 2a; /A A/ 

^dx^- 2T] - HaXMe 



(8) 



By solving Eq. ([5]) with the boundary conditions 
(dri/dx) x= Q — and ^Ix— »oo = 0, we obtain the following 
expression for the depletion of the order parameter: 



r)(x) 



(n\ M e~ 2x/XM - V2e~^ Kx ^ , (9) 



(kX m - V2) 

where rjo is the value of rj(x) at the sample surface, 

H^kXm 



V(0) =Vo = 



2(kX m + V2) ' 



(10) 



This relation is complemented by the expression for the 
magnetic field B z (x), 



B z (x) = H a e- X ' x * 



(11) 



which follows at once from Eq. ©. Both r)(x) and B z {x) 
depend on the variational parameter Am, which can be 
obtained by minimizing the Gibbs free energy [Eq. |T])] . 

In Eq. ©, 770 is related to Am through Eq. ([10]). The 
approximation is more accurate the lower the magnetic 
field is, i.e., when rjo <C 1- In the above equations, we 
have only one variational parameter, which is Am- An 
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alternative possibility is to consider 770 as a second varia- 
tional parameter to obtain a more accurate description of 
the Meissner state up to magnetic fields close to the vor- 
tex penetration field H p . We have followed this second 
procedure in this paper. 

To determine the variational parameters, we must find 
the extremum of given by Eq. |T]), which can be 

written as 



AG 

L Z Ly 



</•' (-/ 2 +^ + i((V/) 2 + u 2 / 2 



/ dx[B(x) - Ha] 
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2 

2 



(12) 



where u is the velocity of the superconducting elec- 
trons, u = (Vip — kA). In the London gauge, u is pro- 
portional to A, u = -kA; therefore, we have 



u x 



0. 



(13) 



The free energy [Eq. (I12| ] can be evaluated by using Eqs. 
©, (HI]), and (H3]) for f{x) = 1 - rj(x), B(x), and u, 
respectively. The minimization of the free energy allows 
us to obtain Am and t]q, which completes the variational 
description of the Meissner state. 



B. Full numerical time-dependent 
Ginzburg— Landau solution 

We compare the variational solutions to the results 
obtained from full numerical solutions of the time- 
dependent GL equations ) 15 ! 16 ' 17 



1 



dt k 
dA _ 1 ^Im[**(V 
~dt ~ V 



2 (V -iKAf-H + (1 -I*! 2 )*, (14) 
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Time is in the units of a characteristic normaliza- 
tion time t — £?/D, where D is the electronic diffu- 
sion constant, a' is the normalized conductivity, cr' = 
(47rA 2 /c 2 r)cr. The other quantities were normalized as in 
the previous section. 

To solve equations (fT4|) and (fl"5|) . we have used the 
standard finite-difference discretization schemed The 
order parameter and vector potential are defined at the 
nodes of a rectangular mesh (r = (/, J)). In our simu- 
lations, we have assumed a sample that is semi-infinite 
in the x direction and infinite in the y, z directions, and 
we have assumed that the magnetic field is applied along 
z. The problem is then reduced to two dimensions be- 
cause we can neglect all derivatives along z. The sym- 
metry of the problem implies that for all mesh points, 
A/, j = (A x i t j, A y i t j, 0) and B/,,/ = {0,Q,B zIiJ ), where 
B zIJ - (V x A), = {d x A yIJ - d y A^j). The link 



variables U j ; j = exp(— inh^A^i.j) (n = x, y) are intro- 
duced in order to preserve gauge invariance in the dis- 
cretization. 

In this geometry, the discretized forms of equations 
(fT4l) and (HSl) are 



0* 


u *i-i,j^i-i,J ~ 2*/, J + 
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(16) 



where Ax and Ay are the mesh widths, cr' was chosen 
as equal to unity, as in Ref. fla . 

The dynamical equations must be complemented with 
the appropriate boundary conditions for both the order 
parameter and the vector potential. We have imposed 
periodic boundary conditions in the y direction, i.e., 

9(x,y) = V(x,y + L v ), 
A x (x,y) = A x {x,y + L y ), 
A y{x,y) = A y {x,y + Ly), 

and semiperiodic boundary conditions in the x direction, 
where one side of the superconductor is in contact with 
the vacuum at x — 0, implying 

((V - zA)*) J^o = 0, 
B\x—q = H a . 

At x = L, we impose the conditions that are obtained 
at x — 00, 



B\ 



By choosing a value much larger than A, L x = 24A 
for L x , we have obtained accurate results for a sample 
semi-infinite in x. 



C. Comparison between the variational solution 
and the full Ginzburg— Landau numerical results 

Figure [1] shows a comparison between the variational 
and full numerical results for the order parameter and 
for the magnetic field in the Meissner state. Both quan- 
tities are calculated along a direction perpendicular to 
the sample surface and for k = 2. The size of the numer- 
ical sample is described by L x — 24 A and L y = 16A. It 
is seen that the variational description is quite accurate, 
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even when H a is near H p , the field of first vortex pene- 
tration [see Figs. Hcl) and[T{c2)]. The numerical simu- 
lations obtain H v = 1.13H C ; this value coincides with the 
results of Ref. 18, which is also a one-dimensional cal- 
culation. In a two-dimensional sample, this result would 
be obtained for a perfect surface that induces a uniform 
depletion of the order parameter along the surface (see 
Ref. [lj^ for a complete discussion) . When a defect^ or 
thermal fluctuations 17 induce the nucleation of a vortex, 
the value of H n diminishes. 
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FIG. 1: Shown is a comparison of the order parameter and 
the magnetic field, in the variational approximation and in the 
numerical calculation, for n — 2 and for different values of the 
applied field. We see that even in (cl) and (c2), when H is 
near the field of first penetration, the variational description 
is quite accurate. 



exact solution by adding the fields produced by the real 
charge and by the image charges. 

Since the GL equations are non-linear, care must be 
exercised in applying the method. To perform a varia- 
tional calculation in the present case, we must physically 
construct acceptable trial functions for the order param- 
eter and currents. 

Clem's^ variational calculation allows us to determine 
the order parameter, field, and current for a vortex in an 
infinite superconductor. In order to use this solution for 
a vortex close to the superconductor-vacuum interface, 
we must first consider a vortex placed at a generic point 
(xo,yo) m a bulk superconductor. We introduce the fol- 
lowing auxiliary variables: 



P (xq, y ; x, y) = \J{x~ x ) 2 + {y- y f 



R (x , Vol x, y) = \J{x- x ) 2 + (y - y ) 2 + Q. 



(19) 



Here, ( v is a variational parameter of the same order 
of magnitude of the coherence length. Clem's variational 
ansatz for the order parameter takes the form 



fvor (x,y) — 



p(x ,y ;x,y) 
R(x ,y ;x,y)' 



(20) 



This allows the exact solution of the second GL equation, 
giving, respectively, for the field and current 



B z = 

u = 



1 K (R(x ,y ;x 1 y)) 

1 p(xo,yo;x,y) Ki(R(x ,y ;x,y)) 



k( v R(x ,y ;x,y) 



Ki(Cv) 



, (21) 



where Kq(x) and Ki(x) are modified Bessel functions. 

As we saw above, the boundary conditions at the sam- 
ple surface require the order parameter to have a zero 
slope there. A second requirement is that no current 
should flow across the sample surface, 



df(x,y) 
dx 



0, 



x=0 

■h (x,y)\ x=0 = 0. 



(22) 



III. CLEM'S VARIATIONAL SOLUTION NEAR 
A SURFACE 

Originally developed for electrostatics and fluid dy- 
namics, the image method was envisaged to automati- 
cally satisfy the boundary conditions in a given problem. 
It has found applications in several fields of physics de- 
scribed by linear equations, wherein the superposition 
principle is valid. In such cases, the method provides the 



Both conditions can be satisfied by considering the com- 
bined effect of the vortex located at the point (xo, yo) plus 
an image vortex, which is located at the point (—xo,yo). 
The currents in the image vortex must rotate in the op- 
posite sense to the ones in the real vortex. The order 
parameter for the image vortex would be 

f , \ P(x,y,-x ,y ) , . 

Jim [X,y,Xo, yo) = — p (23) 

V P (x,y;-xa,y ) 2 + Q 
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To construct the variational order parameter for the 
vortex-image vortex pair, we can simply take the product 

F(x,y;x ,y ) = f V or(x,y;x ,y ) x f im (x, y; -x , y ). 

(24) 

This assumption guarantees the vanishing of the order 
parameter at each vortex core and also of the normal 
slope at the surface, 



dF(x,y) 



dx 



= 0. 



x=0 



When the vortex is placed far enough from the surface, 
Eq. ([24]) tends to the correct limits. 

The variational solution for the current requires some 
care. A velocity field obtained by adding the current 
fields for vortex and image vortex would satisfy the 
boundary condition at the superconductor- vacuum inter- 
face but would violate the requirement that the current 
at each vortex core vanishes. An alternative is to con- 
struct first a compound velocity field by adding the ve- 
locity fields of each vortex. The resulting field would also 
satisfy the boundary condition, with the advantage that 
the singularities at each vortex core would be maintained, 
very much as is the case for the charge-image charge pair 
in electrostatics. The total current field must then be 
calculated from this velocity field. To obtain the veloc- 
ity distribution for the vortex-image vortex system, we 
must consider the sum of the velocities for each of these 
elements as follows: 



U x (x,y)= u x (x,y;x ,y ) -u x (x,y;-x ,yo) , (25) 

U y (x,y)= Uy(x,y;x Q ,y Q )+Uy(x,y;-x ,yo), (26) 
with 

1 K 1 (R(x,y;x ,y )) 



u x {x,y;x ,yo) 



u y {x,y;x ,y ) 



(y-yo)R(x,y;xQ,yg) 

P 2 (x,y;x Q ,y ) 
1 K 1 (R(x,y;xo,y )) 

(x- x ) R(x,y;x ,y ) 
P 2 (x,y;x ,yo) 



In these expressions, u x (x, y; xq, y<j) and 
u y {x,y;—xa,yo) are the velocity field components of a 
vortex centered at (xo,yo), whereas u x [x, y; — Xq, yo) 
and u v (x, y; —xo,yo) are the velocity components of the 
image vortex centered at (— xo,yo). 

This combination keeps the essential property of hav- 
ing the correct divergence at each vortex core, and it 
satisfies the boundary condition 

U x (0,y) = 0. 
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FIG. 2: Variational parameters as a function of the vortex 
position. In (a), parameters rjo and Am, which are related 
to the Meissner state, are shown. Same as in (b) for the 
parameter ( v , which determines the vortex size. System sizes 
are: (•) 20A x 11A, (o) 40A x 11A, and (□) 160A x 11A. 



In order to obtain the variational current field, we must 
combine the velocity field of Eqs. (|25|) and (|26| with the 
order parameter for the vortex-image-vortex pair given 
by Eq. (pH)) . The total current would thus be 



J(x,y) = \F(x,y)\ 2 V(x,y). 

It can be seen that this expression satisfies the correct 
boundary condition as follows: 

J*(0,y) = 0. 

The total magnetic field of the vortex-antivortex pair 
would be 

B z (x, y) = B z (x, y; x ,y ) - B z (x, y; -x , y ) 

or 

n / s _ J_ f K (R(x,y;x ,y )) _ K (R (x,y; -x ,y )) 

z[x,y) ~< v { ^(C) ifi(C«) 

(27) 

To evaluate the free energy given by Eq. (fT2)l . we 
have to combine the vortex-antivortex expressions with 
the contributions due to the Meissner currents. For the 
magnetic field and currents, we assume a superposition 
principle and simply add the contributions due to each 
source. Thus, to obtain the total magnetic field, we have 
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to add Eqs. ([TTJ) and (|27|). The velocity of the super- 
conducting electrons u is obtained by adding the vector 
component given by Eq. (fT3| and Eqs. (|25p and (|26|) . 
For the order parameter, instead, we have to multiply 
both contributions, given by Eq. (|24[) and by the contri- 
bution of the Meissner state, f(x) = 1 — 77(2), with r){x) 
given by Eq. The minimization of the free energy 

allows us to obtain £„, Am , and 770, which completes the 
variational description of a vortex near a surface. 



A. Steadiness of the variational parameters 

To test the steadiness of the variational calculation, 
we have studied the change of the variational parameters 
as a function of £0, the vortex position. We have also 
checked the convergence of these parameters in terms of 
the size of the numerical sample. Figure [5] shows the 
behavior of the three variational parameters, the vortex 
size £„, the Meissner parameter tjq, and Am, as functions 
of the vortex position for different sample sizes. 

It can be seen in Fig. &[a,) that with an increase in 
the sample size, the last two parameters become inde- 
pendent of the vortex position, which is an indication of 
the adequacy of the variational function. On the other 
hand, the vortex parameter ( v is strongly affected when 
the vortex moves close to the surface [Fig. (Hb)] and 
does not show appreciable changes with increasing sam- 
ple size. The value of for large x coincides with the 
value obtained by Clem for k = 2 in Ref. [H, ( v — 1.15£. 

The results of Fig. [2] show that we can obtain the 
Meissner parameters 770 and Am from the Meissner vari- 
ational calculation and use them as fixed parameters in 
the vortex variational equations. This allows a faster 
convergence of the solution in the presence of vortices 
because we only need to minimize the energy with re- 
spect to a single parameter, the vortex core size £«• We 
have used this procedure in what follows. 

Our results show good agreement between the varia- 
tional calculations and full numerical results both for the 
energy barrier, as will be shown in the next section, and 
for other quantities, particularly the shape of the vortex 
near the surface. 



B. Quality of the variational solution of a vortex 
near a surface 

As we did in Sec. II for the Meissner state, in this 
section we compare the variational solutions for a vortex 
near a surface to the results obtained from the full nu- 
merical solutions of the GL equations. However, in the 
present case, the comparison needs to be more carefully 
done. In the variational calculation, the position of the 
vortex at a given point, xq, is fixed and the energy must 
be minimized in order to find the variational parameters. 
Due to the forces exerted by the Meissner currents, such 
a configuration is unstable in the GL case. The diffi- 




FIG. 3: (Color online) Comparison between the profiles of the 
order parameter from (a) the variational model and (b) full 
numerical solution of the GL equations. Both graphs are for 
k = 2, x = 2.9A, and H a = 0.70H C . 

culties in pinning an isolated vortex at position xq are 
not solely related to the time dependent equations that 
we are using to find the equilibrium configurations. In 
a time independent approach, the system also tends to 
the equilibrium configuration that for a field larger than 
the field of first vortex penetration (H p ) is a vortex at 
position xq — > 00. 

In order to overcome this difficulty and to calculate the 
energy of a vortex located at Xq, we pin the vortex by 
using a square numerical seed of size d = £ for the order 
parameter. The use of a pinning seed poses the problem 
of the distortion of the order parameter around the vortex 
core, which affects the evaluation of the free energy. We 
reduce this side effect by introducing a seed that has 
the same shape as a vortex in bulk. The solutions of the 
numerical equations then converge towards a stable state 
containing a vortex pinned at position Xq. We move the 
seed in steps given by the spatial discretization, allowing 
the system to relax to a new stable solution at each new 
position of the seed. Our choice of pinning seed allows a 
good comparison to the variational solution of a vortex 
at position xq. A similar procedure was used in Ref. 0, 
wherein by fixing the phase of the order parameter, it 
was possible to pin and move vortices near a surface. 

Figure [3] shows a comparison between the profiles of the 
order parameter obtained from the variational method 
and from the full numerical simulations of the GL equa- 
tions. The figures are for xo = 2.9A and H a = Q.70H C . 
As we see from Figs. EJa) anddjb), the description ob- 
tained from the variational model is quite accurate, al- 
though some qualitative differences are apparent. This 
is a typical feature of variational calculations, wherein 
generally a better agreement is obtained for the energy 
calculation than for the field properties. 

In Fig. |4l we show a comparison between the varia- 
tional calculations and numerical results for the energy 
barrier as a function of the vortex position and for dif- 
ferent values of the applied magnetic field. The maxima 
of the energy as a function of the vortex position for 
fields lower than H p generate the energy barrier for vor- 
tex penetration. The energy barrier can be defined as 
the energy difference between the value at the maximum 



7 



=2; H a =0.42H c I (a) 



-0.409- 
-0.410- 
j -0.41 1 
-0.412- 
-0.413- 
-0.414- 
-0.415- 



ff \ H — ■ — Numerical 
Pi ■□ — p— Variational 

/ ^ 




t \ 

ff ±\ — ■— Numerical 

I \ D \ — o— Variational 



x IX 



(b) 



- 



(c) 



x IX 



X 



x IX 



x IX 



'**«( 




x/X 



FIG. 4: Shown is a comparison between variational calcula- 
tion and numerical results for the energy barrier as a function 
of the vortex position for k — 2 and for different values of 
the applied field. Arrows follow the increase in magnetic field 
from (a) to (d). It is seen that even in (d), wherein H a is near 
Hp = l.OiHc, the variational description is quite accurate. 



and the value at the surface, A = G(x max ) — G(0). By 
increasing the magnetic field from H a = 0A2H C [Fig. 
Ufa)] to H a = 0.70F C [Fig. gfc)], it is seen that the 
energy barrier decreases and, finally, almost disappears 
near H a = 0.85H C , as shown in Fig. HJd). At the same 
time, the maximum moves closer to the surface. 

It is seen that even when H a is near H pi the variational 
description is quite accurate [see Figs. 0Jc) and 2] (d)]. 
These results should be compared to those shown in Fig. 
3 of Ref. Il3j, wherein it is seen that the energy of the 
vortex line quite accurately coincides with the full GL 
numerical results in a wide range of k values. 

For k = 2, the field of first penetration obtained 
turns out to be H p = 1.0AH C , which is a value lower 
than the one obtained in Sec. II for the Meissner state 
(Hp = 1.13H C ). In the full numerical simulations, this 
is a consequence of the symmetry breaking produced by 
the pinning center we have used; while in the variational 
approach, the symmetry is already broken by the nature 
of the solution we have imposed. 

In Figs. HJc) andUJd), the energy differences between 
both approaches are higher near the surface, xq w 0. 
This is a consequence of the increase in the force that the 



FIG. 5: Same as Fig. gjbut for k = 3. 

Meissner currents exert on the vortex. In the full numer- 
ical approach, this means a higher difficulty in pinning 
a vortex at position close to xq ~ 0. In any case, the 
energy differences are always lower than 0.05%. 

Figure O shows a similar comparison to that in Fig. 
[4] in the case of n — 3. As can be seen from Figs. [5jc) 
and[5]^d), when H a is near H p , the variational description 
remains as accurate as in the previous case. The max- 
ima of the energy as a function of the vortex position for 
fields lower than H p generates the energy barrier for vor- 
tex penetration. By increasing the magnetic field from 
H a = 0A2H C [Fig. EKa)] to H a = 0.68ff c [Fig. EJc)], the 
energy barrier decreases and, finally, almost disappears 
near H a = Q.8H C , as shown in Fig. [5jd) . For n = 3, the 
field of first penetration turns out to be H p = 0.91H C . 

The usefulness of the variational approach can be 
stressed by calculating other quantities related to vor- 
tices near surfaces. One such quantity is the magnetic 
flux. As early as 1961, Bardeer>2i showed that magnetic 
flux in a superconducting cylinder can be less than one 
flux quantum. Later in Ref. I22I , Shmidt and Mkrtchyan, 
using an extension of the London model, calculated the 
magnetic flux for a vortex near the surface of a semi- 
infinite sample. They found the following functional de- 
pendence: 

$ = I B-da = $o(l - e~ X0 ). (28) 
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1 2 3 4 5 6 

X Q /l 

FIG. 6: Magnetic flux of a vortex as a function of the distance 
to the surface. The continuous curve corresponds to Eq. l|28[) 
and the squares to the variational calculation. 

We have used the variational model to calculate the 
magnetic flux as a function of the distance to the sample 
surface, as shown in Fig. [6] Our results agree quite well 
with the functional dependence of Eq. (|28p . We note 
that the fluxoid quantization 

e»> 

remains valid, even when $ is less than one flux quantum, 
due to the contribution of the superconducting currents 

Js- 

Geim et al., using a Hall probe in Ref. d experimentally 
confirmed the fact that vortices can have less than one 
flux quantum in mesoscopic samples. Similar results had 
previously been obtained in experiments on bulk samples 
by Civale and de la Crua 2 ^. In Ref. [H, they studied 
the magnetic behavior as a function of temperature of 
samples with a constant number of vortices pinned at 
fixed positions. They observed that the magnetic flux 
carried by vortices located close to the surface increases 
with decreasing temperature due to an indirect increase 
in the distance to the surface when A(T) decreases. 



IV. EXTENSION TO TWO VORTICES 

The variational model of Clem was previously ex- 
tended to describe a flux lattice in Ref. [24|, wherein the 
reversible magnetization of high-T c superconductors as a 
function of the applied field was calculated. Our formula- 
tion can be straightforwardly extended to the description 
of two or more vortices near a surface. In particular, we 
focus on the case of two vortices, located at a distance Xo 
from the sample surface and separated by a distance yo- 
In order to calculate the interaction force, we first fol- 
low the same procedure used in Sec. Ill to calculate the 



variational energy of the system. We introduce an order 
parameter, which is the product of Clem's variational 
expressions for the two vortices and the corresponding 
images. Similarly, the velocities and the total magnetic 
field are obtained by following the procedure described 
in Sec. III. Once the energy of the system containing two 
vortices is obtained, the interaction force between them 
can be calculated from the numerical derivative of the 
energy of the system. 

We concentrate first on the interaction force between 
vortices that are away from the sample surface, i.e. for 
xq — > oo. In this case, the image vortices can be omitted 
from the calculations because their influence is negligi- 
ble when they are far from the surface. In Fig. we 
show the interaction force for two cases: Fig. UJa) is 
for k — 10 and Fig. Etb) is for k = 2. A comparison is 
shown with the force calculated within the London model 
and also with the formula obtained from the long-range 
asymptotic behavio r 25 i 26 within the Ginzburg-Landau 
approach. Good agreement between all curves is ob- 
tained for distances larger than 0.5A for k = 10 and for 
distances larger than 2. OA for k = 2. It should be noted 
that in the case of the London model, the force corre- 
sponding to two vortices Fl = (An/ K 2 )K (ri — rj) di- 
verges when the distance between vortices (r^ — Tj ) tends 
to zero, which in our case is when yo — > 0. 

The expression for the long-range asymptotic behavior 
of the GL model, due to Kramer, i o 25 ' 26 

F int = -J {tfofa - r 3 ) - K [V2K{n - rj )]}. (30) 

This expression incorporates a correction due to the 
overlapping of the vortex cores that induces an attractive 
component in the force between vortices. The final result 
gives a finite value for the force corresponding to the long- 
range asymptotic GL model when yo — > 0, as can be seen 
in Figs. Eta) and[?tb). 

On the other hand, the attractive contribution in the 
variational model results in a zero force for yo = when 
the two vortices merge in a two-quanta vortex. This con- 
figuration is unstable in a type II superconductor, and a 
minimal separation between vortices would result in a re- 
pulsive force. An interaction force decreasing and going 
to zero for very small vortex separation is in agreement 
with previous variational calculations of the interaction 
energy between vortices obtained by Jacobs and Rebbi 
in Ref. ,27. As we can see in Fig. the interaction force 
obtained from the variational model has a maximum lo- 
cated at yo = 0.4A for k = 10 and at yo = 1-3A when 
k = 2. 

In Fig. [5J we fix the distance yo between two vor- 
tices and calculate the variation of the interaction force 
between them as a function of the distance xo to the sur- 
face. In this case, the contribution of the image vortices 
becomes more important as the pair approaches the sur- 
face. We show the interaction force parallel to the surface 
(along the line that connects both vortices) as a function 
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FIG. 7: Force between two vortices as a function of the dis- 
tance between them when both vortices are well apart from 
the surface. In (a) and (b), we show a comparison of the be- 
havior in the London (A), long-range GL (v)> an< A variational 
(■) models. 



of xq for k = 2 in Fig. Eta) and for k = 10 in Fig. Etb). 
In Fig. Eta), the open squares are for H a = 0.4 and the 
closed squares are for H a = 0; in both cases, we observe 
a steady decrease in the force between vortices when the 
distance to the surface decreases. There is only a small 
difference between the forces obtained at different values 
of the applied magnetic field due to the differences be- 
tween the Meissner currents induced in each case. The 
overall qualitative behavior is independent of the value 
of the distance between vortices yo, as we see in Fig. Eta) 
by comparing the results for j/o = 2 A and j/o — 2.5 A. In 
Fig. EJb), the same qualitative behavior is observed for 
k = 10. From Fig. El we can conclude that there is a 
steady decrease in the interaction force between vortices 
when the distance to the surface decreases and that the 
force goes to zero when the pair is close to the boundary. 

This result can be understood by considering that 
the repulsive interaction force between two vortices is 
screened by the contribution of the attractive interaction 
force between each vortex with the image of the other 
vortex. In particular, when xq — > 0, a vortex and its im- 
age are located at approximately the same position for 



the other vortex and the interaction force goes to zero. 

Our results suggest that the vortex lattice is softer in 
the direction parallel to the surface in finite samples. It 
should be interesting to devise experiments that can ex- 
plore these properties. Even when we have calculated 
the case of a semi-infinite sample, the same qualitative 
behavior is expected to appear in a thin film with vor- 
tices parallel to the surface. In a thin film, vortices are 
confined by two surfaces and the image vortices corre- 
sponding to both surfaces contribute to the screening of 
the vortex interaction forces. 
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FIG. 8: Interaction force between two vortices separated a 
distance j/o as a function of their distance to the surface xq. 
In (a) and (b), we compare results obtained at different values 
of re, yo, and applied magnetic field H a - 



V. CONCLUSIONS 

We have shown that Clem's variational ansatz for a 
free vortex can be extended to the description of vortex 
penetration. The results show quite good agreement with 
the full numerical results both for the energy barrier and 
for the description of the vortex near the surface. The 
flux carried by a vortex as a function of its distance to 
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the surface can be shown to be easily calculated and to 
coincide with known results. 

We extended the model to calculate the force between 
two vortices. When the vortices are far from the sur- 
face, the variational results show good agreement with 
the London and long-range GL results for large intervor- 
tex distances; whereas for small distances, the variational 
model gives vanishing forces corresponding to the merg- 
ing of the two vortices in a double quantized vortex^. 
We also found a steady decrease in the interaction force 
between vortices when the distance to the surface de- 
creases; the interaction force goes to zero when the pair 
is close to the boundary. 

Our variational approach gives manageable expressions 
that can be used to obtain approximations for all phys- 
ically relevant quantities. Another advantage of this 
method is the lower computational time that it requires, 
allowing one to obtain fast and reliable solutions for a 
vortex near a surface. The agreement between the varia- 
tional solution and numerical calculations shows the use- 
fulness of the former for intermediate k when computa- 
tions become heavy. For large k, numerical calculations 
can be based on the London model for the magnetic con- 
tribution. This description is not accurate at lower k 
values where a numerical approach to the GL descrip- 
tion is more appropriate. A particular problem arises at 
intermediate k, at which the computation within the GL 
model becomes very demanding because of the difficulty 
to describe both spatial scales with the same discretiza- 
tion. It is in this range that the variational approach is 
most welcome. 

The method used in this paper can be generalized for 



mesoscopic superconductors. However, consideration of 
more than one surface in some cases could give rise to 
infinite images. This difficulty can be overcome by trun- 
cating the infinite series as was done in Ref. [29j for vor- 
tex penetration in a thin film using the London model. 
In this paper, we have assumed a semi-infinite medium 
with no demagnetizing effects. In this case, the boundary 
condition for the magnetic field B| s — H a applies at the 
sample surface. A similar condition applies for a thin film 
with the externally applied magnetic field parallel to the 
surface. A thin film is an example of a mesoscopic sys- 
tem wherein our results can be generalized in a straight- 
forward manner. However, in general, demagnetization 
effects are important in mesoscopic superconductors of 
finite thickness where the boundary condition B| s = H a 
applies at infinity and not at the sample boundary. Vor- 
tices in mesoscopic superconductors are confined by the 
sample surface and their interaction with the Meissner 
currents is very important, a situation with similarities 
to the case analyzed in this work. This opens up inter- 
esting questions about the behavior of the effective inter- 
action force between two vortex cores as a function of the 
distance to the surface in mesoscopic superconductors. 
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